!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the subroutine is used to calcaulate point on symmetry axis  
! papramenters by MOC.
! ref. to chapter 16 of 'Gas Dynamics' by Maurice J. Zucrow & Joe D. Hoffmann
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-06-18
! revised  :
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! include  :
!   subroutine : SymmetryP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   1. angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine SymmetryP()
   use VariableDef
   implicit none
!     1
!           4()
!     (2)
!        given variable
!   real*8::x1,y1,Vu1,Vv1
   integer::icor,i
!        temp variable
   real*8::Vum,Vvm,ym,Vm,thetam,am,Mm,lm,Qm,Rm,Sm,Tm
   real*8::x4f,Vu4f
!        calculate variable
!   real*8::x4,Vu4,V4,p4,T4,rou4,M4
   
   icor=6
   x4f=0.0
   Vu4f=0.0
   y1=abs(y1)
   Vv1=abs(Vv1)
   Vum=Vu1
   Vvm=Vv1
   ym=y1

   do i=1,icor
      Vm=sqrt(Vum**2+Vvm**2)
      thetam=atan(Vvm/Vum)
      call Thermo(Vm,p4,rou4,T4,am,Mm)
      
      lm=tan(thetam-asin(1/Mm))
      Qm=Vum**2-am**2
      Rm=2*Vum*Vvm-Qm*lm
      Sm=delta*Vvm*am**2/ym
      
      x4=x1-y1/lm
      y4=0.0   ! x-axis in the symmetry
      
      Tm=Sm*(x4-x1)+Qm*Vu1+Rm*Vv1
      
      Vu4=Tm/Qm
      Vv4=0.0   ! axial-symmetry point

      if (abs(x4-x4f) <= el) then
         if (abs(Vu4-Vu4f) <= ev) then
            exit
         end if
      end if

      Vum=(Vu1+Vu4)/2
      Vvm=Vv1/2
      ym=y1/2
      
      x4f=x4
      Vu4f=Vu4
   end do
   V4=sqrt(Vu4**2+Vv4**2)
   call Thermo(V4,p4,rou4,T4,am,M4)
   
   if (wallflag==-1) then
      y1=-y1
      Vv1=-Vv1
   end if
end subroutine SymmetryP
